####-------------------------------------------------------------------------------------####
# This file analyzes  2008 panel

####-------------------------------------------------------------------------------------####
rm(list=ls())
detach("package:car")
require(foreign)
require(car)
###  Load Data ##
load("/Users/chrisweber/Dropbox/Working Projects/Nick Davis_Polarization Emotion/data/panel2008FINAL.RData")
d1<-panel.2008

##Emotions
d1$anxiety.out.w1<-(rowMeans(cbind(d1$anger.out.w1, d1$fear.out.w1), na.rm=T)-1)/4
d1$enthusiasm.in.w1<-(rowMeans(cbind(d1$hope.in.w1, d1$pride.in.w1), na.rm=T)-1)/4
d1$anxiety.in.w1<-(rowMeans(cbind(d1$anger.in.w1, d1$fear.in.w1), na.rm=T)-1)/4
d1$enthusiasm.out.w1<-(rowMeans(cbind(d1$hope.out.w1, d1$proud.out.w1), na.rm=T)-1)/4

d1$anxiety.out.w11<-(rowMeans(cbind(d1$anger.out.w11, d1$fear.out.w11), na.rm=T)-1)/4
d1$enthusiasm.in.w11<-(rowMeans(cbind(d1$hope.in.w11, d1$proud.in.w11), na.rm=T)-1)/4
d1$anxiety.in.w11<-(rowMeans(cbind(d1$anger.in.w11, d1$fear.in.w11), na.rm=T)-1)/4
d1$enthusiasm.out.w11<-(rowMeans(cbind(d1$hope.out.w11, d1$proud.out.w11), na.rm=T)-1)/4
table(d1$pid.w1)
table(d1$ideo.w1)
## Sorting Measure
d1$sorting.w1<-recode(abs(d1$pid.w1-d1$ideo.w1)+1, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1")*(abs(4-d1$pid.w1)+0)*(abs(4-d1$ideo.w1)+0)
d1$sorting.w11<-recode(abs(d1$pid.w11-d1$ideo.w11)+1, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1")*(abs(4-d1$pid.w11)+0)*(abs(4-d1$ideo.w11)+0)

table(d1$sorting.w1)
table(d1$sorting.w11)
####################################################################################
####################################################################################
##

# Main Sorting Models: 
m1a<-lm(sorting.w11~
          sorting.w1+anxiety.in.w1+anxiety.out.w1+enthusiasm.in.w1+enthusiasm.out.w1+
          catholic.2008+other.2008+
          sex.2008+college.2008+news+nonwhite.2008+
          age.2008, d1) 
m1b<-lm(anxiety.in.w11~
          sorting.w1+anxiety.in.w1+anxiety.out.w1+enthusiasm.in.w1+enthusiasm.out.w1+
          catholic.2008+other.2008+
          sex.2008+college.2008+news+nonwhite.2008+
          age.2008, d1) 
m1c<-lm(anxiety.out.w11~
          sorting.w1+anxiety.in.w1+anxiety.out.w1+enthusiasm.in.w1+enthusiasm.out.w1+
          catholic.2008+other.2008+
          sex.2008+college.2008+news+nonwhite.2008+
          age.2008, d1) 
m1d<-lm(enthusiasm.in.w11~
          sorting.w1+anxiety.in.w1+anxiety.out.w1+enthusiasm.in.w1+enthusiasm.out.w1+
          catholic.2008+other.2008+
          sex.2008+college.2008+news+nonwhite.2008+
          age.2008, d1) 
m1e<-lm(enthusiasm.out.w11~
          sorting.w1+anxiety.in.w1+anxiety.out.w1+enthusiasm.in.w1+enthusiasm.out.w1+
          catholic.2008+other.2008+
          sex.2008+college.2008+news+nonwhite.2008+
          age.2008, d1) 


output.coef1<-summary(m1a)$coefficients[,2]
output.se1<-summary(m1a)$coefficients[,3]

output.coef2<-summary(m1b)$coefficients[,2]
output.se2<-summary(m1b)$coefficients[,3]

output.coef3<-summary(m1c)$coefficients[,2]
output.se3<-summary(m1c)$coefficients[,3]

output.coef4<-summary(m1d)$coefficients[,2]
output.se4<-summary(m1d)$coefficients[,3]

output.coef5<-summary(m1e)$coefficients[,2]
output.se5<-summary(m1e)$coefficients[,3]

require(dplyr)
output<-cbind(output.coef1, output.se1,
              output.coef2, output.se2,
              output.coef3, output.se3,
              output.coef4, output.se4,
              output.coef5, output.se5)%>%
  as.data.frame()
names(output)<-rep(c("Slope", "t"), times=5)

require(xtable)
xtable(output, digits=3,
       caption="")


summary(m1a)$r.squared
summary(m1b)$r.squared
summary(m1c)$r.squared
summary(m1d)$r.squared
summary(m1e)$r.squared

dim(model.matrix(m1a))[1]
dim(model.matrix(m1b))[1]
dim(model.matrix(m1c))[1]
dim(model.matrix(m1d))[1]
dim(model.matrix(m1e))[1]


##### Estimate the Markov Model ######
d2<-panel.2000
d2$id<-c(1:dim(d2)[1])

m1e<-lm(enthusiasm.out.w11~
          sorting.w1+anxiety.in.w1+anxiety.out.w1+enthusiasm.in.w1+enthusiasm.out.w1+
          catholic.2008+other.2008+
          sex.2008+college.2008+news+nonwhite.2008+
          age.2008, d1) 



d2$anxiety.out<-(rowMeans(cbind(d2$anger.out, d2$fear.out), na.rm=T)-1)/4
d2$enthusiasm.in<-(rowMeans(cbind(d2$hope.in, d2$pride.in), na.rm=T)-1)/4
d2$anxiety.in<-(rowMeans(cbind(d2$anger.in, d2$fear.in), na.rm=T)-1)/4
d2$enthusiasm.out<-(rowMeans(cbind(d2$hope.out, d2$proud.out), na.rm=T)-1)/4

d2$ideology<-(d2$ideo.2000-1)/6
d2$intAnxIn<-d2$ideology*d2$anxiety.in
d2$intAnxOut<-d2$ideology*d2$anxiety.out
d2$intEnthIn<-d2$ideology*d2$enthusiasm.in
d2$intEnthOut<-d2$ideology*d2$enthusiasm.out



data<-data.frame(id=rep(1:length(d2$id), times=3), state=c(d2$pid.3.2000, d2$pid.3.2002, d2$pid.3.2004),
                 ideology=rep(d2$ideology, times=3),
                 sex=rep(d2$sex.2000, times=3),
                 college=rep(d2$college.2000, times=3),
                 age=(rep(d2$age.2000, times=3)-18)/79,
                 catholic=rep(d2$catholic.2000, times=3),
                 other=rep(d2$other.2000, times=3),
                 anxiety.in=rep(d2$anxiety.in, times=3),
                 anxiety.out=rep(d2$anxiety.out, times=3),
                 enthusiasm.in=rep(d2$enthusiasm.in, times=3),
                 enthusiasm.out=rep(d2$enthusiasm.out, times=3),
                 intAnxIn=rep(d2$intAnxIn, times=3),
                 intAnxOut=rep(d2$intAnxOut, times=3),
                 intEnthIn=rep(d2$intEnthIn, times=3),
                 intEnthOut=rep(d2$intEnthOut, times=3),
                time=c(rep(c(1:3), each=length(d2$id) )))
require(dplyr)
data<-data%>%
  arrange(id, time)%>%
  na.omit
statetable.msm(state, id, data)
twoway3.q<-rbind(c(0,0.10, 0),
                 c(0.1,0, 0.10),
                 c(0, 0.10, 0)
)
#rownames(twoway3.q)<-colnames(twoway3.q)<-c("Democrat", "Independent", "Republican")
twoway3.q
model<-msm(state~time, covariates=~ideology+anxiety.in+anxiety.out+enthusiasm.in+enthusiasm.out+intAnxIn
           +intAnxOut+intEnthIn+intEnthOut, 
           subject=id, data=data,
           qmatrix=twoway3.q, gen.inits=TRUE, obstype=1,
           method="BFGS", control=list(trace=1, REPORT=1, fnscale=2000, maxit=20000, reltol=1e-8))

hazard.msm(model)

### Consistent Affect, Liberals

t1<-pmatrix.msm(model, t=2, ci="none", covariates=list(ideology=0, 
                                                        anxiety.in=0,
                                                        anxiety.out=1,
                                                        enthusiasm.in=0,
                                                        enthusiasm.out=0,
                                                        intAnxOut=0,
                                                        intAnxIn=0,
                                                        intEnthIn=0,
                                                        intEnthOut=0
))


print(t1)
print(o)




d1$liberal.s<-ifelse(d1$id.3cat=="Liberal", 1, 0)
d1$conservative.s<-ifelse(d1$id.3cat=="Conservative", 1, 0)
d1$moderate.s<-ifelse(d1$id.3cat=="Moderate", 1, 0)

d1$intAnxIn1<-d1$liberal.s*d1$anxiety.in
d1$intAnxOut1<-d1$liberal.s*d1$anxiety.out
d1$intEnthIn1<-d1$liberal.s*d1$enthusiasm.in
d1$intEnthOu1t<-d1$liberal.s*d1$enthusiasm.out

d1$intAnxIn2<-d1$conservative.s*d1$anxiety.in
d1$intAnxOut2<-d1$conservative.s*d1$anxiety.out
d1$intEnthIn2<-d1$conservative.s*d1$enthusiasm.in
d1$intEnthOut2<-d1$conservative.s*d1$enthusiasm.out

 # Sam's heatmap
 {
# # Sam's Heatmap
# map<-polr(as.factor(pid)~
#             ideology3+anxiety.in+anxiety.out+enthusiasm.in+enthusiasm.out+
#             intAnxIn+intAnxOut+intEnthIn+intEnthOut+
#             catholic+no.religion+other.religion+racial.resentment+
#             female+education+news+non.white+
#             age+year2012+year2016, data=d1) 
# ## Vary Anxiety, hold constant everthing else
# ## Vary Enthusiasm, hold constant 
# emotion1<-seq(0, 1, by=0.1)
# emotion2<-seq(0,1, by=0.1)
# covariate.grid<-expand.grid(emotion1=emotion1, emotion2=emotion2)
# 
# contour.fun<-function(output, emotion1=0, emotion2=0, ideology="liberal", target="consistent"){
#   require(MASS)
#   t<-model.matrix(output)[,-1]
#   if(ideology=="liberal" & target=="consistent"){
#     p<-cbind(
#       0,
#       t[,2],
#       emotion1,
#       emotion2,
#       t[,5],
#       0*t[,2],
#       0*emotion1,
#       0*emotion2,
#       0*t[,5],
#       t[,10:20]
#     )
# }
#   if(ideology=="conservative" & target=="consistent"){
#     p<-cbind(
#       1,
#       t[,2],
#       emotion1,
#       emotion2,
#       t[,5],
#       1*t[,2],
#       1*emotion1,
#       1*emotion2,
#       1*t[,5],
#       t[,10:20]
#     )
#   }
#   if(ideology=="liberal" & target=="inconsistent"){
#     p<-cbind(
#       0,
#       emotion1,
#       t[,3],
#       t[,4],
#       emotion2,
#       0*emotion1,
#       0*t[,3],
#       0*t[,4],
#       0*emotion2,
#       t[,10:20]
#     )
#   }
#   if(ideology=="conservative" & target=="inconsistent"){
#     p<-cbind(
#       1,
#       emotion1,
#       t[,3],
#       t[,4],
#       emotion2,
#       1*emotion1,
#       1*t[,3],
#       1*t[,4],
#       1*emotion2,
#       t[,10:20]
#     )
#   }
#     g1<-apply(cbind(p, -1, 0, 0, 0, 0, 0)%*%c(output$coefficients, output$zeta), 2, mean) ### The average effect. Simulate b and auth
#     g2<-apply(cbind(p, 0, -1, 0, 0, 0, 0)%*%c(output$coefficients, output$zeta),2, mean) ### The average effect. Simulate b and auth
#     g3<-apply(cbind(p, 0, 0, -1, 0, 0, 0)%*%c(output$coefficients, output$zeta),2, mean) ### The average effect. Simulate b and auth
#     g4<-apply(cbind(p, 0, 0, 0,-1, 0, 0)%*%c(output$coefficients, output$zeta),2, mean) ### The average effect. Simulate b and auth
#     g5<-apply(cbind(p, 0, 0, 0, 0, -1, 0)%*%c(output$coefficients, output$zeta),2, mean) ### The average effect. Simulate b and auth
#     g6<-apply(cbind(p, 0, 0, 0, 0, 0, -1)%*%c(output$coefficients, output$zeta),2, mean) ### The average effect. Simulate b and auth
#   
#     pr.1<-(1-plogis(g3))###Democrat
#     pr.2<-(plogis(g3)-plogis(g4))  ##Independent
#     pr.3<-plogis(g4) ###Republcan
#     
#     return(list(Democrat=data.frame(pr.1),
#                 Independent=data.frame(pr.2),
#                 Republican= data.frame(pr.3)
#     ))
# }
# 
# prediction1<-contour.fun(map, covariate.grid[1,1], covariate.grid[1,2], "liberal", "consistent")[[1]]
# prediction2<-contour.fun(map, covariate.grid[1,1], covariate.grid[1,2], "liberal", "inconsistent")[[1]]
# prediction3<-contour.fun(map, covariate.grid[1,1], covariate.grid[1,2], "conservative", "consistent")[[3]]
# prediction4<-contour.fun(map, covariate.grid[1,1], covariate.grid[1,2], "conservative", "inconsistent")[[3]]
# 
# 
# for(i in 2:dim(covariate.grid)[1]){
#   prediction1[i]<-rbind(contour.fun(map, covariate.grid[i,1], covariate.grid[i,2], "liberal", "consistent"))[[1]]
#   prediction2[i]<-rbind(contour.fun(map, covariate.grid[i,1], covariate.grid[i,2], "liberal", "inconsistent"))[[1]]
#   prediction3[i]<-rbind(contour.fun(map, covariate.grid[i,1], covariate.grid[i,2], "conservative", "consistent"))[[3]]
#   prediction4[i]<-rbind(contour.fun(map, covariate.grid[i,1], covariate.grid[i,2], "conservative", "inconsistent"))[[3]]
# }
# plot.data<-data.frame(cbind(liberal1=as.numeric(prediction1), 
#                             liberal2=as.numeric(prediction2), 
#                             conservative1=as.numeric(prediction3), 
#                             conservative2=as.numeric(prediction4), 
#                             negative=covariate.grid[,1], 
#                             positive=covariate.grid[,2]))
# 
# 
# require(akima)
# 
# pdf("/users/chrisweber/Dropbox/Project Folder/projects/Nick Davis_Polarization Emotion/paper/sam1.pdf", width=8, height=7)
# par(mfcol=c(1,2)) ## create outer margin
# z1<-interp(plot.data$positive,  plot.data$negative, plot.data$liberal1)
# image(z1,  col=gray((20:0)/20),   xlab="Positive Affect", ylab="Negative Affect", plot.title=title(main="Liberal"), 
#       useRaster=TRUE, cex.axis=0.6, cex.lab=0.6, cex.main=0.3)
# contour(z1, add=T, col = "white", nlevels=15)
# 
# z1<-interp(plot.data$positive,  plot.data$negative, plot.data$conservative1)
# image(z1,  col=gray((20:0)/20),   xlab="Positive Affect", ylab="Negative Affect", plot.title=title(main="Conservative"), 
#       useRaster=TRUE, cex.axis=0.6, cex.lab=0.6, cex.main=0.3)
# contour(z1, add=T, col = "white", nlevels=15)
# mtext("Consistent Emotions",  outer = TRUE, cex=1.2)
# dev.off
# par(mfcol=c(1,2)) ## create outer margin
# 
# pdf("/users/chrisweber/Dropbox/Project Folder/projects/Nick Davis_Polarization Emotion/paper/sam2.pdf", width=8, height=7)
# z1<-interp(plot.data$positive,  plot.data$negative, plot.data$liberal2)
# image(z1,  col=gray((20:0)/20),   xlab="Positive Affect", ylab="Negative Affect", plot.title=title(main="Liberal"), 
#       useRaster=TRUE, cex.axis=0.6, cex.lab=0.6, cex.main=0.3)
# contour(z1, add=T, col = "white", nlevels=15)
# 
# z1<-interp(plot.data$positive,  plot.data$negative, plot.data$conservative2)
# image(z1,  col=gray((20:0)/20),   xlab="Positive Affect", ylab="Negative Affect", plot.title=title(main="Conservative"), 
#       useRaster=TRUE, cex.axis=0.6, cex.lab=0.6, cex.main=0.3)
# contour(z1, add=T, col = "white", nlevels=15)
# mtext("Inconsistent Emotions",  outer = TRUE, cex=1.2)
# dev.off()
# 
# 
# 
# 
# z1<-interp(plot.data$positive,  plot.data$negative, plot.data$liberal1)
# image(z1,  col=gray((20:0)/20),   xlab="Negative Affect", ylab="Positive Affect", plot.title=title(main="Pr(Dem)|Liberal"))
# contour(z1, add=T, col = "white", nlevels=15)
# 
# 
# 
# z1<-interp(plot.data$positive,  plot.data$negative, plot.data$liberal1)
# image(z1,  col=gray((20:0)/20),   xlab="Comparative Anxiety", ylab="Ideology (Conservative)", plot.title=title(main="Pr(Identify as Republican)"))
# contour(z1, add=T, col = "white", nlevels=15)
# 
# z1<-interp(plot.data$positive,  plot.data$negative, plot.data$liberal1)
# image(z1,  col=gray((20:0)/20),   xlab="Comparative Anxiety", ylab="Ideology (Conservative)", plot.title=title(main="Pr(Identify as Republican)"))
# contour(z1, add=T, col = "white", nlevels=15)
# 
# 
# z1<-interp(plot.data$emotion,  plot.data$ideology, log(plot.data$enthusiasm))
# image(z1,  col=gray((30:0)/30), xlab="Comparative Enthusiasm", ylab="Ideology (Conservative)", plot.title=title(main="Pr(Identify as Republican)"))
# contour(z1, add=T, col = "white", nlevels=15)
# 
 }
 

m1a<-polr(as.factor(pid)~
           ideology3+anxiety.in+anxiety.out+enthusiasm.in+enthusiasm.out+
            intAnxIn+
           catholic+no.religion+other.religion+racial.resentment+
           female+education+news+non.white+
           age+year2012+year2016, data=d1)  
m1b<-polr(as.factor(pid)~
            ideology3+anxiety.in+anxiety.out+enthusiasm.in+enthusiasm.out+
            intAnxOut+
            catholic+no.religion+other.religion+racial.resentment+
            female+education+news+non.white+
            age+year2012+year2016, data=d1)  
m1c<-polr(as.factor(pid)~
            ideology3+anxiety.in+anxiety.out+enthusiasm.in+enthusiasm.out+
            intEnthIn+
            catholic+no.religion+other.religion+racial.resentment+
            female+education+news+non.white+
            age+year2012+year2016, data=d1)  
m1d<-polr(as.factor(pid)~
            ideology3+anxiety.in+anxiety.out+enthusiasm.in+enthusiasm.out+
            intEnthIn+
            catholic+no.religion+other.religion+racial.resentment+
            female+education+news+non.white+
            age+year2012+year2016, data=d1)  
m2<-polr(as.factor(pid)~
           ideology3+anxiety.in+anxiety.out+enthusiasm.in+enthusiasm.out+intAnxIn+
           intAnxOut+intEnthIn+intEnthOut+
           catholic+no.religion+other.religion+racial.resentment+
           female+education+news+non.white+
           age+year2012+year2016, data=d1)  


  
# PID/ Table 1/ Figure 2
## Table ###
out1<-rbind(coef(summary(m1a)) [c(1:6),1:2], c(NA,NA), c(NA,NA), c(NA,NA)) #ANES
out2<-rbind(coef(summary(m1b)) [c(1:5),1:2], c(NA,NA), coef(summary(m1b)) [c(6),1:2], c(NA,NA), c(NA,NA)) #ANES
out3<-rbind(coef(summary(m1c)) [c(1:5),1:2], c(NA,NA), c(NA,NA), coef(summary(m1c)) [c(6),1:2], c(NA,NA)) #ANES
out4<-rbind(coef(summary(m1d)) [c(1:5),1:2], c(NA,NA), c(NA,NA), c(NA,NA), coef(summary(m1d)) [c(6),1:2]) #ANES
out5<-rbind(coef(summary(m2)) [c(1:9),1:2]) #ANES


results1<-as.matrix(
  cbind(out1, out2, out3, out4, out5
  ))

row.names(results1)<-c("Ideology", "In Party Negative Affect", "Out Party Negative Affect", "In Party Positive Affect", "Out Party Positive Affect", "Out Party Negative x Ideology", "In Party Negative x Ideology", "Out Party Positive x Ideology", "In Party Positive x Ideology")
require(xtable)                  
xtable(results1, digits=3,
       caption="Partisan-Ideological Sorting by Emotion. Ordinal logistic regression estimates. Entries are point estimates and standard errors. 
       Entries in bold indicate a slope that is two times the estimated standard error.
       All models control for religion, racial resentment, gender, education, cross section, news consumption, race and age. The full models
       may be found in the appendix.")


xtable(summary(m2)$coefficients, digits=3,
       caption="Partisan-Ideological Sorting by Emotion. Ordinal logistic regression estimates. Entries are point estimates and standard errors. 
       Entries in bold indicate a slope that is two times the estimated standard error.
       All models control for religion, racial resentment, gender, education, cross section, news consumption, race and age. The full models
       may be found in the appendix.")



#########################################################################################################
#########################################################################################################
##################### Figure 1 ###############################################################
#########################################################################################################
## Functions to plot data ##
## Figure plots three plots -- marginal effects for ideology across four groups x pid

# pred.generator is the function that generates the predicted marginal effects.
marginal.effect.plot<-function(party="Democrat"){
  seq.number<-0.1
  value<-seq(0,1, seq.number)
  if(party=="Democrat"){
enth.in<-pred.generator(m2, "in.party.enthusiasm", 0, "ologit")[[1]]
enth.out<-pred.generator(m2, "out.party.enthusiasm", 0, "ologit")[[1]]
anx.in<-pred.generator(m2, "in.party.anxiety", 0, "ologit")[[1]]
anx.out<-pred.generator(m2, "out.party.anxiety", 0, "ologit")[[1]]

for(i in 2:length(value)){
  enth.in[i,]<-rbind(pred.generator(m2, "in.party.enthusiasm", value[i], "ologit")[[1]])
  enth.out[i,]<-rbind(pred.generator(m2, "out.party.enthusiasm", value[i], "ologit")[[1]])
  anx.in[i,]<-rbind(pred.generator(m2, "in.party.anxiety", value[i], "ologit")[[1]])
  anx.out[i,]<-rbind(pred.generator(m2, "out.party.anxiety", value[i], "ologit")[[1]])
}

enth.in$emotion<-seq(0, 1, seq.number)
enth.in$type<-"In Party Enthusiasm"
enth.out$emotion<-seq(0, 1, seq.number)
enth.out$type<-"Out Party Enthusiasm"
anx.in$emotion<-seq(0, 1, seq.number)
anx.in$type<-"In Party Anxiety"
anx.out$emotion<-seq(0, 1, seq.number)
anx.out$type<-"Out Party Anxiety"
}

  if(party=="Independent"){
    enth.in<-pred.generator(m2, "in.party.enthusiasm", value=0, "ologit")[[2]]
    enth.out<-pred.generator(m2, "out.party.enthusiasm", 0, "ologit")[[2]]
    anx.in<-pred.generator(m2, "in.party.anxiety", 0, "ologit")[[2]]
    anx.out<-pred.generator(m2, "out.party.anxiety", 0, "ologit")[[2]]
    
    for(i in 2:length(value)){
      enth.in[i,]<-rbind(pred.generator(m2, "in.party.enthusiasm", value[i], "ologit")[[2]])
      enth.out[i,]<-rbind(pred.generator(m2, "out.party.enthusiasm", value[i], "ologit")[[2]])
      anx.in[i,]<-rbind(pred.generator(m2, "in.party.anxiety", value[i], "ologit")[[2]])
      anx.out[i,]<-rbind(pred.generator(m2, "out.party.anxiety", value[i], "ologit")[[2]])
    }
    
    enth.in$emotion<-seq(0, 1, seq.number)
    enth.in$type<-"In Party Enthusiasm"
    enth.out$emotion<-seq(0, 1, seq.number)
    enth.out$type<-"Out Party Enthusiasm"
    anx.in$emotion<-seq(0, 1, seq.number)
    anx.in$type<-"In Party Anxiety"
    anx.out$emotion<-seq(0, 1, seq.number)
    anx.out$type<-"Out Party Anxiety"
  }
  
  if(party=="Republican"){
    enth.in<-pred.generator(m2, "in.party.enthusiasm", 0, "ologit")[[3]]
    enth.out<-pred.generator(m2, "out.party.enthusiasm", 0, "ologit")[[3]]
    anx.in<-pred.generator(m2, "in.party.anxiety", 0, "ologit")[[3]]
    anx.out<-pred.generator(m2, "out.party.anxiety", 0, "ologit")[[3]]
    
    for(i in 2:length(value)){
      enth.in[i,]<-rbind(pred.generator(m2, "in.party.enthusiasm", value[i], "ologit")[[3]])
      enth.out[i,]<-rbind(pred.generator(m2, "out.party.enthusiasm", value[i], "ologit")[[3]])
      anx.in[i,]<-rbind(pred.generator(m2, "in.party.anxiety", value[i], "ologit")[[3]])
      anx.out[i,]<-rbind(pred.generator(m2, "out.party.anxiety", value[i], "ologit")[[3]])
    }
    
    enth.in$emotion<-seq(0, 1, seq.number)
    enth.in$type<-"In Party Enthusiasm"
    enth.out$emotion<-seq(0, 1, seq.number)
    enth.out$type<-"Out Party Enthusiasm"
    anx.in$emotion<-seq(0, 1, seq.number)
    anx.in$type<-"In Party Anxiety"
    anx.out$emotion<-seq(0, 1, seq.number)
    anx.out$type<-"Out Party Anxiety"
  }
  return(rbind(enth.in, enth.out, anx.in, anx.out))
}
plot.data.democrat<-marginal.effect.plot("Democrat")
plot.data.republican<-marginal.effect.plot("Republican")
plot.data.independent<-marginal.effect.plot("Independent")

in.party.enthusiasm<-rbind(subset(plot.data.democrat, type=="In Party Enthusiasm"),
      subset(plot.data.republican, type=="In Party Enthusiasm"),
      subset(plot.data.independent, type=="In Party Enthusiasm")
      )
in.party.enthusiasm$Party<-rep(c("Democrat", "Republican", "Independent"), each=11)

out.party.enthusiasm<-rbind(subset(plot.data.democrat, type=="Out Party Enthusiasm"),
                           subset(plot.data.republican, type=="Out Party Enthusiasm"),
                           subset(plot.data.independent, type=="Out Party Enthusiasm")
)

out.party.enthusiasm$Party<-rep(c("Democrat", "Republican", "Independent"), each=11)

in.party.anxiety<-rbind(subset(plot.data.democrat, type=="In Party Anxiety"),
                           subset(plot.data.republican, type=="In Party Anxiety"),
                           subset(plot.data.independent, type=="In Party Anxiety")
)
in.party.anxiety$Party<-rep(c("Democrat", "Republican", "Independent"), each=11)

out.party.anxiety<-rbind(subset(plot.data.democrat, type=="Out Party Anxiety"),
                            subset(plot.data.republican, type=="Out Party Anxiety"),
                            subset(plot.data.independent, type=="Out Party Anxiety")
)
out.party.anxiety$Party<-rep(c("Democrat", "Republican", "Independent"), each=11)




### Generate the Figure 
names(in.party.enthusiasm)<-c("min1", "min2", "max2", "max1", "mean", "emotion", "type", "Party")
names(out.party.anxiety)<-c("min1", "min2", "max2", "max1", "mean", "emotion", "type", "Party")
names(in.party.anxiety)<-c("min1", "min2", "max2", "max1", "mean", "emotion", "type", "Party")
names(out.party.enthusiasm)<-c("min1", "min2", "max2", "max1", "mean", "emotion", "type", "Party")

{
plot1<-
  ggplot(out.party.anxiety) + 
  geom_ribbon(aes(x=emotion, ymin=min1, ymax=max1, group=Party), fill="grey", alpha=0.5)+
  geom_ribbon(aes(x=emotion, ymin=min2, ymax=max2, group=Party), fill="grey", alpha=0.5)+
  scale_colour_manual(name="Party", values=c("blue", "black", "red"))+
  geom_line(aes(x=emotion, y=mean, colour=Party))+
  theme_bw() +
  #  theme(legend.position="none") +
  theme(panel.background=element_rect(fill="#F0F0F0")) +
  theme(plot.background=element_rect(fill="#F0F0F0")) +
  theme(panel.border=element_rect(colour="#F0F0F0")) +
  # Format the grid
  theme(panel.grid.major=element_line(colour="#D0D0D0",size=.25)) +
  theme(axis.ticks=element_blank())+
  ggtitle("Out Party Negative Affect") +
  theme(plot.title=element_text(face="bold",hjust=-.08,vjust=2,colour="#3C3C3C",size=11)) +
    theme(axis.text.x=element_text(size=9,colour="#535353",face="bold")) +
    theme(axis.text.y=element_text(size=9,colour="#535353",face="bold")) +
    theme(axis.title.y=element_text(size=9,colour="#535353",face="bold",vjust=1.5)) +
    theme(axis.title.x=element_text(size=9,colour="#535353",face="bold",vjust=-.5)) +
    scale_y_continuous("Marginal Effect of Ideology", limits=c(-1,1))+
    scale_x_continuous("Emotion")+
  geom_hline(yintercept=0)


plot2<-
  ggplot(in.party.enthusiasm) + 
  geom_ribbon(aes(x=emotion, ymin=min1, ymax=max1, group=Party), fill="grey", alpha=0.5)+
  geom_ribbon(aes(x=emotion, ymin=min2, ymax=max2, group=Party), fill="grey", alpha=0.5)+
  scale_colour_manual(name="Party", values=c("blue", "black", "red"))+
  geom_line(aes(x=emotion, y=mean, colour=Party))+
  theme_bw() +
  #  theme(legend.position="none") +
  theme(panel.background=element_rect(fill="#F0F0F0")) +
  theme(plot.background=element_rect(fill="#F0F0F0")) +
  theme(panel.border=element_rect(colour="#F0F0F0")) +
  # Format the grid
  theme(panel.grid.major=element_line(colour="#D0D0D0",size=.25)) +
  theme(axis.ticks=element_blank())+
  ggtitle("In Party Positive Affect") +
  theme(plot.title=element_text(face="bold",hjust=-.08,vjust=2,colour="#3C3C3C",size=11)) +
  theme(axis.text.x=element_text(size=9,colour="#535353",face="bold")) +
  theme(axis.text.y=element_text(size=9,colour="#535353",face="bold")) +
  theme(axis.title.y=element_text(size=9,colour="#535353",face="bold",vjust=1.5)) +
  theme(axis.title.x=element_text(size=9,colour="#535353",face="bold",vjust=-.5)) +
  scale_y_continuous("Marginal Effect of Ideology", limits=c(-1,1))+
  scale_x_continuous("Emotion")+
  geom_hline(yintercept=0)


plot3<-
  ggplot(in.party.anxiety) + 
  geom_ribbon(aes(x=emotion, ymin=min1, ymax=max1, group=Party), fill="grey", alpha=0.5)+
  geom_ribbon(aes(x=emotion, ymin=min2, ymax=max2, group=Party), fill="grey", alpha=0.5)+
  scale_colour_manual(name="Party", values=c("blue", "black", "red"))+
  geom_line(aes(x=emotion, y=mean, colour=Party))+
  theme_bw() +
  #  theme(legend.position="none") +
  theme(panel.background=element_rect(fill="#F0F0F0")) +
  theme(plot.background=element_rect(fill="#F0F0F0")) +
  theme(panel.border=element_rect(colour="#F0F0F0")) +
  # Format the grid
  theme(panel.grid.major=element_line(colour="#D0D0D0",size=.25)) +
  theme(axis.ticks=element_blank())+
  ggtitle("In Party Negative Affect") +
  theme(plot.title=element_text(face="bold",hjust=-.08,vjust=2,colour="#3C3C3C",size=11)) +
  theme(axis.text.x=element_text(size=9,colour="#535353",face="bold")) +
  theme(axis.text.y=element_text(size=9,colour="#535353",face="bold")) +
  theme(axis.title.y=element_text(size=9,colour="#535353",face="bold",vjust=1.5)) +
  theme(axis.title.x=element_text(size=9,colour="#535353",face="bold",vjust=-.5)) +
  scale_y_continuous("Marginal Effect of Ideology", limits=c(-1,1))+
  scale_x_continuous("Emotion")+
  geom_hline(yintercept=0)


plot4<-
  ggplot(out.party.enthusiasm) + 
  geom_ribbon(aes(x=emotion, ymin=min1, ymax=max1, group=Party), fill="grey", alpha=0.5)+
  geom_ribbon(aes(x=emotion, ymin=min2, ymax=max2, group=Party), fill="grey", alpha=0.5)+
  scale_colour_manual(name="Party", values=c("blue", "black", "red"))+
  geom_line(aes(x=emotion, y=mean, colour=Party))+
  theme_bw() +
  #  theme(legend.position="none") +
  theme(panel.background=element_rect(fill="#F0F0F0")) +
  theme(plot.background=element_rect(fill="#F0F0F0")) +
  theme(panel.border=element_rect(colour="#F0F0F0")) +
  # Format the grid
  theme(panel.grid.major=element_line(colour="#D0D0D0",size=.25)) +
  theme(axis.ticks=element_blank())+
  ggtitle("Out Party Positive Affect") +
  theme(plot.title=element_text(face="bold",hjust=-.08,vjust=2,colour="#3C3C3C",size=11)) +
  theme(axis.text.x=element_text(size=9,colour="#535353",face="bold")) +
  theme(axis.text.y=element_text(size=9,colour="#535353",face="bold")) +
  theme(axis.title.y=element_text(size=9,colour="#535353",face="bold",vjust=1.5)) +
  theme(axis.title.x=element_text(size=9,colour="#535353",face="bold",vjust=-.5)) +
  scale_y_continuous("Marginal Effect of Ideology", limits=c(-1,1))+
  scale_x_continuous("Emotion")+
  geom_hline(yintercept=0)
}

  pdf("/users/chrisweber/Google Drive/Project Folder/projects/Nick Davis_Polarization Emotion/paper/graph4.pdf", width=8, height=7)
  grid_arrange_shared_legend(plot1, plot2, plot3, plot4, ncol=2, nrow=2)
  dev.off()
  ####################################################################################
  #################################################################################
  ##################################################  ##########################################
##### Marginal Effects for Consistency######
  ##################################################  ##########################################
  ####################################################################################
  # PID/ Table 1/ Figure 2
  ## Table ###
  d1$SintAnxIn3<-d1$social.ideology*d1$anxiety.in
  d1$SintAnxOut3<-d1$social.ideology*d1$anxiety.out
  d1$SintEnthIn3<-d1$social.ideology*d1$enthusiasm.in
  d1$SintEnthOut3<-d1$social.ideology*d1$enthusiasm.out

  d1$FintAnxIn3<-d1$fiscal.ideology*d1$anxiety.in
  d1$FintAnxOut3<-d1$fiscal.ideology*d1$anxiety.out
  d1$FintEnthIn3<-d1$fiscal.ideology*d1$enthusiasm.in
  d1$FintEnthOut3<-d1$fiscal.ideology*d1$enthusiasm.out
  

  m1<-polr(as.factor(pid)~social.ideology+fiscal.ideology+anxiety.in+anxiety.out+
             enthusiasm.in+enthusiasm.out+SintAnxIn3+
             SintAnxOut3+SintEnthIn3+SintEnthOut3+
             FintAnxIn3+
             FintAnxOut3+FintEnthIn3+FintEnthOut3+
             catholic+no.religion+other.religion+racial.resentment+
             female+education+news+non.white+
             age+year2012+year2016, data=d1) ## Same here ##
  
  m2<-lm(social.ideology~fiscal.ideology+anxiety.in+anxiety.out+
             enthusiasm.in+enthusiasm.out+FintAnxIn3+
             FintAnxOut3+FintEnthIn3+FintEnthOut3+
             catholic+no.religion+other.religion+racial.resentment+
             female+education+news+non.white+
             age+year2012+year2016, data=d1) ## Same here ##
  
  
  xtable(summary(m1)$coefficients[,1:3], digits=3,
         caption="Ideological Consistency. Logistic regression estimates. 
          Entries are point estimates, standard errors, and t-values.")
  
  xtable(summary(m2)$coefficients, digits=3,
         caption="Ideological Consistency. Ordered logistic regression estimates. 
          Entries are point estimates, standard errors, and t-values.")
  

  anx.in<-predicted.value.logit(m1, 0, emotion="anxiety.in", model="logit")
  anx.out<-predicted.value.logit(m1, 0, emotion="anxiety.out", model="logit")
  enth.in<-predicted.value.logit(m1, 0, emotion="enthusiasm.in", model="logit")
  enth.out<-predicted.value.logit(m1, 0, emotion="enthusiasm.out", model="logit")
  value<-seq(0,1, 0.1)
  for(i in 2:length(value)){
    enth.in[i,]<-rbind(predicted.value.logit(m1, value[i], emotion="enthusiasm.in", model="logit"))
    enth.out[i,]<-rbind(predicted.value.logit(m1, value[i],emotion="enthusiasm.out",  model="logit"))
    anx.in[i,]<-rbind(predicted.value.logit(m1, value[i],emotion="anxiety.in", model="logit"))
    anx.out[i,]<-rbind(predicted.value.logit(m1,value[i], emotion="anxiety.out",  model="logit"))
  }    

  names(enth.in)<-c("min1", "min2", "max2", "max1", "mean")
  names(enth.out)<-c("min1", "min2", "max2", "max1", "mean")
  names(anx.in)<-c("min1", "min2", "max2", "max1", "mean")
  names(anx.out)<-c("min1", "min2", "max2", "max1", "mean")
  enth.in$emotion<-value
  enth.out$emotion<-value
  anx.in$emotion<-value
  anx.out$emotion<-value
  
  

  {
    plot1<-
      ggplot(enth.in) + 
      geom_ribbon(aes(x=emotion, ymin=min1, ymax=max1), fill="grey", alpha=0.5)+
      geom_line(aes(x=emotion, y=mean))+
      theme_bw() +
      #  theme(legend.position="none") +
      theme(panel.background=element_rect(fill="#F0F0F0")) +
      theme(plot.background=element_rect(fill="#F0F0F0")) +
      theme(panel.border=element_rect(colour="#F0F0F0")) +
      # Format the grid
      theme(panel.grid.major=element_line(colour="#D0D0D0",size=.25)) +
      theme(axis.ticks=element_blank())+
      ggtitle("In Party Positive Affect") +
      theme(plot.title=element_text(face="bold",hjust=-.08,vjust=2,colour="#3C3C3C",size=10)) +
      theme(axis.text.x=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.text.y=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.title.y=element_text(size=9,colour="#535353",face="bold",vjust=1.5)) +
      theme(axis.title.x=element_text(size=9,colour="#535353",face="bold",vjust=-.5)) +
      scale_y_continuous("Pr(Consistent)", limits=c(0,0.5))+
      scale_x_continuous("Emotion")

    plot2<-
      ggplot(enth.out) + 
      geom_ribbon(aes(x=emotion, ymin=min1, ymax=max1), fill="grey", alpha=0.5)+
      geom_line(aes(x=emotion, y=mean))+
      theme_bw() +
      #  theme(legend.position="none") +
      theme(panel.background=element_rect(fill="#F0F0F0")) +
      theme(plot.background=element_rect(fill="#F0F0F0")) +
      theme(panel.border=element_rect(colour="#F0F0F0")) +
      # Format the grid
      theme(panel.grid.major=element_line(colour="#D0D0D0",size=.25)) +
      theme(axis.ticks=element_blank())+
      ggtitle("Out Party Positive Affect") +
      theme(plot.title=element_text(face="bold",hjust=-.08,vjust=2,colour="#3C3C3C",size=10)) +
      theme(axis.text.x=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.text.y=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.title.y=element_text(size=9,colour="#535353",face="bold",vjust=1.5)) +
      theme(axis.title.x=element_text(size=9,colour="#535353",face="bold",vjust=-.5)) +
      scale_y_continuous("Pr(Consistent)", limits=c(0,0.5))+
      scale_x_continuous("Emotion")
    
    
    plot3<-
      ggplot(anx.in) + 
      geom_ribbon(aes(x=emotion, ymin=min1, ymax=max1), fill="grey", alpha=0.5)+
      geom_line(aes(x=emotion, y=mean))+
      theme_bw() +
      #  theme(legend.position="none") +
      theme(panel.background=element_rect(fill="#F0F0F0")) +
      theme(plot.background=element_rect(fill="#F0F0F0")) +
      theme(panel.border=element_rect(colour="#F0F0F0")) +
      # Format the grid
      theme(panel.grid.major=element_line(colour="#D0D0D0",size=.25)) +
      theme(axis.ticks=element_blank())+
      ggtitle("In Party Negative Affect") +
      theme(plot.title=element_text(face="bold",hjust=-.08,vjust=2,colour="#3C3C3C",size=10)) +
      theme(axis.text.x=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.text.y=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.title.y=element_text(size=9,colour="#535353",face="bold",vjust=1.5)) +
      theme(axis.title.x=element_text(size=9,colour="#535353",face="bold",vjust=-.5)) +
      scale_y_continuous("Pr(Consistent)", limits=c(0,0.5))+
      scale_x_continuous("Emotion")
    
    plot4<-
      ggplot(anx.out) + 
      geom_ribbon(aes(x=emotion, ymin=min1, ymax=max1), fill="grey", alpha=0.5)+
      geom_line(aes(x=emotion, y=mean))+
      theme_bw() +
      #  theme(legend.position="none") +
      theme(panel.background=element_rect(fill="#F0F0F0")) +
      theme(plot.background=element_rect(fill="#F0F0F0")) +
      theme(panel.border=element_rect(colour="#F0F0F0")) +
      # Format the grid
      theme(panel.grid.major=element_line(colour="#D0D0D0",size=.25)) +
      theme(axis.ticks=element_blank())+
      ggtitle("Out Party Negative Affect") +
      theme(plot.title=element_text(face="bold",hjust=-.08,vjust=2,colour="#3C3C3C",size=10)) +
      theme(axis.text.x=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.text.y=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.title.y=element_text(size=9,colour="#535353",face="bold",vjust=1.5)) +
      theme(axis.title.x=element_text(size=9,colour="#535353",face="bold",vjust=-.5)) +
      scale_y_continuous("Pr(Consistent)", limits=c(0,0.5))+
      scale_x_continuous("Emotion")
    
  }
  
  pdf("/users/chrisweber/Google Drive/Project Folder/projects/Nick Davis_Polarization Emotion/paper/graph_consistent.pdf", width=8, height=7)
  multiplot(plot4, plot2, plot1, plot3, cols=2)
  dev.off()
  
  
marginal.effect.ideology<-function(party="Liberal"){
    seq.number<-0.1
    value<-seq(0,1, seq.number)
    if(party=="Liberal"){
      enth.in<-pred.generator(m2, "in.party.enthusiasm", 0, "ologit")[[1]]
      enth.out<-pred.generator(m2, "out.party.enthusiasm", 0, "ologit")[[1]]
      anx.in<-pred.generator(m2, "in.party.anxiety", 0, "ologit")[[1]]
      anx.out<-pred.generator(m2, "out.party.anxiety", 0, "ologit")[[1]]
      
      for(i in 2:length(value)){
        enth.in[i,]<-rbind(pred.generator(m2, "in.party.enthusiasm", value[i], "ologit")[[1]])
        enth.out[i,]<-rbind(pred.generator(m2, "out.party.enthusiasm", value[i], "ologit")[[1]])
        anx.in[i,]<-rbind(pred.generator(m2, "in.party.anxiety", value[i], "ologit")[[1]])
        anx.out[i,]<-rbind(pred.generator(m2, "out.party.anxiety", value[i], "ologit")[[1]])
      }
      
      enth.in$emotion<-seq(0, 1, seq.number)
      enth.in$type<-"In Party Enthusiasm"
      enth.out$emotion<-seq(0, 1, seq.number)
      enth.out$type<-"Out Party Enthusiasm"
      anx.in$emotion<-seq(0, 1, seq.number)
      anx.in$type<-"In Party Anxiety"
      anx.out$emotion<-seq(0, 1, seq.number)
      anx.out$type<-"Out Party Anxiety"
    }
    
    if(party=="Moderate"){
      enth.in<-pred.generator(m2, "in.party.enthusiasm", 0, "ologit")[[2]]
      enth.out<-pred.generator(m2, "out.party.enthusiasm", 0, "ologit")[[2]]
      anx.in<-pred.generator(m2, "in.party.anxiety", 0, "ologit")[[2]]
      anx.out<-pred.generator(m2, "out.party.anxiety", 0, "ologit")[[2]]
      
      for(i in 2:length(value)){
        enth.in[i,]<-rbind(pred.generator(m2, "in.party.enthusiasm", value[i], "ologit")[[2]])
        enth.out[i,]<-rbind(pred.generator(m2, "out.party.enthusiasm", value[i], "ologit")[[2]])
        anx.in[i,]<-rbind(pred.generator(m2, "in.party.anxiety", value[i], "ologit")[[2]])
        anx.out[i,]<-rbind(pred.generator(m2, "out.party.anxiety", value[i], "ologit")[[2]])
      }
      
      enth.in$emotion<-seq(0, 1, seq.number)
      enth.in$type<-"In Party Enthusiasm"
      enth.out$emotion<-seq(0, 1, seq.number)
      enth.out$type<-"Out Party Enthusiasm"
      anx.in$emotion<-seq(0, 1, seq.number)
      anx.in$type<-"In Party Anxiety"
      anx.out$emotion<-seq(0, 1, seq.number)
      anx.out$type<-"Out Party Anxiety"
    }
    
    if(party=="Conservative"){
      enth.in<-pred.generator(m2, "in.party.enthusiasm", 0, "ologit")[[3]]
      enth.out<-pred.generator(m2, "out.party.enthusiasm", 0, "ologit")[[3]]
      anx.in<-pred.generator(m2, "in.party.anxiety", 0, "ologit")[[3]]
      anx.out<-pred.generator(m2, "out.party.anxiety", 0, "ologit")[[3]]
      
      for(i in 2:length(value)){
        enth.in[i,]<-rbind(pred.generator(m2, "in.party.enthusiasm", value[i], "ologit")[[3]])
        enth.out[i,]<-rbind(pred.generator(m2, "out.party.enthusiasm", value[i], "ologit")[[3]])
        anx.in[i,]<-rbind(pred.generator(m2, "in.party.anxiety", value[i], "ologit")[[3]])
        anx.out[i,]<-rbind(pred.generator(m2, "out.party.anxiety", value[i], "ologit")[[3]])
      }
      
      enth.in$emotion<-seq(0, 1, seq.number)
      enth.in$type<-"In Party Enthusiasm"
      enth.out$emotion<-seq(0, 1, seq.number)
      enth.out$type<-"Out Party Enthusiasm"
      anx.in$emotion<-seq(0, 1, seq.number)
      anx.in$type<-"In Party Anxiety"
      anx.out$emotion<-seq(0, 1, seq.number)
      anx.out$type<-"Out Party Anxiety"
    }
    return(rbind(enth.in, enth.out, anx.in, anx.out))
  }
  plot.data.liberal<-marginal.effect.ideology("Liberal")
  plot.data.conservative<-marginal.effect.ideology("Conservative")
  plot.data.moderate<-marginal.effect.ideology("Moderate")
  
  in.party.enthusiasm<-rbind(subset(plot.data.liberal, type=="In Party Enthusiasm"),
                             subset(plot.data.moderate, type=="In Party Enthusiasm"),
                             subset(plot.data.conservative, type=="In Party Enthusiasm")
  )
  in.party.enthusiasm$Party<-rep(c("Liberal", "Moderate", "Conservative"), each=11)
  
  out.party.enthusiasm<-rbind(subset(plot.data.liberal, type=="Out Party Enthusiasm"),
                              subset(plot.data.moderate, type=="Out Party Enthusiasm"),
                              subset(plot.data.conservative, type=="Out Party Enthusiasm")
  )
  
  out.party.enthusiasm$Party<-rep(c("Liberal", "Moderate", "Conservative"), each=11)
  
  in.party.anxiety<-rbind(subset(plot.data.liberal, type=="In Party Anxiety"),
                          subset(plot.data.moderate, type=="In Party Anxiety"),
                          subset(plot.data.conservative, type=="In Party Anxiety")
  )
  in.party.anxiety$Party<-rep(c("Liberal", "Moderate", "Conservative"), each=11)
  
  out.party.anxiety<-rbind(subset(plot.data.liberal, type=="Out Party Anxiety"),
                           subset(plot.data.moderate, type=="Out Party Anxiety"),
                           subset(plot.data.conservative, type=="Out Party Anxiety")
  )
  out.party.anxiety$Party<-rep(c("Liberal", "Moderate", "Conservative"), each=11)
  
  ### Generate the Figure 
  names(in.party.enthusiasm)<-c("min1", "min2", "max2", "max1", "mean", "emotion", "type", "Party")
  names(out.party.anxiety)<-c("min1", "min2", "max2", "max1", "mean", "emotion", "type", "Party")
  names(in.party.anxiety)<-c("min1", "min2", "max2", "max1", "mean", "emotion", "type", "Party")
  names(out.party.enthusiasm)<-c("min1", "min2", "max2", "max1", "mean", "emotion", "type", "Party")
  
  {
    plot1<-
      ggplot(out.party.anxiety) + 
      geom_ribbon(aes(x=emotion, ymin=min1, ymax=max1, group=Party), fill="grey", alpha=0.5)+
      geom_ribbon(aes(x=emotion, ymin=min2, ymax=max2, group=Party), fill="grey", alpha=0.5)+
      scale_colour_manual(name="Ideology", values=c("red", "blue", "black"))+
      geom_line(aes(x=emotion, y=mean, colour=Party))+
      theme_bw() +
      #  theme(legend.position="none") +
      theme(panel.background=element_rect(fill="#F0F0F0")) +
      theme(plot.background=element_rect(fill="#F0F0F0")) +
      theme(panel.border=element_rect(colour="#F0F0F0")) +
      # Format the grid
      theme(panel.grid.major=element_line(colour="#D0D0D0",size=.25)) +
      theme(axis.ticks=element_blank())+
      ggtitle("Out Party Negative Affect") +
      theme(plot.title=element_text(face="bold",hjust=-.08,vjust=2,colour="#3C3C3C",size=11)) +
      theme(axis.text.x=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.text.y=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.title.y=element_text(size=9,colour="#535353",face="bold",vjust=1.5)) +
      theme(axis.title.x=element_text(size=9,colour="#535353",face="bold",vjust=-.5)) +
      scale_y_continuous("Marginal Effect of Operational Ideology", limits=c(-1,1))+
      scale_x_continuous("Emotion")+
      geom_hline(yintercept=0)
    
    
    plot2<-
      ggplot(in.party.enthusiasm) + 
      geom_ribbon(aes(x=emotion, ymin=min1, ymax=max1, group=Party), fill="grey", alpha=0.5)+
      geom_ribbon(aes(x=emotion, ymin=min2, ymax=max2, group=Party), fill="grey", alpha=0.5)+
      scale_colour_manual(name="Ideology", values=c("red", "blue", "black"))+
      geom_line(aes(x=emotion, y=mean, colour=Party))+
      theme_bw() +
      #  theme(legend.position="none") +
      theme(panel.background=element_rect(fill="#F0F0F0")) +
      theme(plot.background=element_rect(fill="#F0F0F0")) +
      theme(panel.border=element_rect(colour="#F0F0F0")) +
      # Format the grid
      theme(panel.grid.major=element_line(colour="#D0D0D0",size=.25)) +
      theme(axis.ticks=element_blank())+
      ggtitle("In Party Positive Affect") +
      theme(plot.title=element_text(face="bold",hjust=-.08,vjust=2,colour="#3C3C3C",size=11)) +
      theme(axis.text.x=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.text.y=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.title.y=element_text(size=9,colour="#535353",face="bold",vjust=1.5)) +
      theme(axis.title.x=element_text(size=9,colour="#535353",face="bold",vjust=-.5)) +
      scale_y_continuous("Marginal Effect of Operational Ideology", limits=c(-1,1))+
      scale_x_continuous("Emotion")+
      geom_hline(yintercept=0)
    
    
    plot3<-
      ggplot(in.party.anxiety) + 
      geom_ribbon(aes(x=emotion, ymin=min1, ymax=max1, group=Party), fill="grey", alpha=0.5)+
      geom_ribbon(aes(x=emotion, ymin=min2, ymax=max2, group=Party), fill="grey", alpha=0.5)+
      scale_colour_manual(name="Ideology", values=c("red", "blue", "black"))+
      geom_line(aes(x=emotion, y=mean, colour=Party))+
      theme_bw() +
      #  theme(legend.position="none") +
      theme(panel.background=element_rect(fill="#F0F0F0")) +
      theme(plot.background=element_rect(fill="#F0F0F0")) +
      theme(panel.border=element_rect(colour="#F0F0F0")) +
      # Format the grid
      theme(panel.grid.major=element_line(colour="#D0D0D0",size=.25)) +
      theme(axis.ticks=element_blank())+
      ggtitle("In Party Negative Affect") +
      theme(plot.title=element_text(face="bold",hjust=-.08,vjust=2,colour="#3C3C3C",size=11)) +
      theme(axis.text.x=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.text.y=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.title.y=element_text(size=9,colour="#535353",face="bold",vjust=1.5)) +
      theme(axis.title.x=element_text(size=9,colour="#535353",face="bold",vjust=-.5)) +
      scale_y_continuous("Marginal Effect of Operational Ideology", limits=c(-1,1))+
      scale_x_continuous("Emotion")+
      geom_hline(yintercept=0)
    
    
    plot4<-
      ggplot(out.party.enthusiasm) + 
      geom_ribbon(aes(x=emotion, ymin=min1, ymax=max1, group=Party), fill="grey", alpha=0.5)+
      geom_ribbon(aes(x=emotion, ymin=min2, ymax=max2, group=Party), fill="grey", alpha=0.5)+
      scale_colour_manual(name="Ideology", values=c("red", "blue", "black"))+
      geom_line(aes(x=emotion, y=mean, colour=Party))+
      theme_bw() +
      #  theme(legend.position="none") +
      theme(panel.background=element_rect(fill="#F0F0F0")) +
      theme(plot.background=element_rect(fill="#F0F0F0")) +
      theme(panel.border=element_rect(colour="#F0F0F0")) +
      # Format the grid
      theme(panel.grid.major=element_line(colour="#D0D0D0",size=.25)) +
      theme(axis.ticks=element_blank())+
      ggtitle("Out Party Positive Affect") +
      theme(plot.title=element_text(face="bold",hjust=-.08,vjust=2,colour="#3C3C3C",size=11)) +
      theme(axis.text.x=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.text.y=element_text(size=9,colour="#535353",face="bold")) +
      theme(axis.title.y=element_text(size=9,colour="#535353",face="bold",vjust=1.5)) +
      theme(axis.title.x=element_text(size=9,colour="#535353",face="bold",vjust=-.5)) +
      scale_y_continuous("Marginal Effect of Operational Ideology", limits=c(-1,1))+
      scale_x_continuous("Emotion")+
      geom_hline(yintercept=0)
  }
  
  pdf("/users/chrisweber/Google Drive/Project Folder/projects/Nick Davis_Polarization Emotion/paper/graph_ideology.pdf", width=8, height=7)
  grid_arrange_shared_legend(plot1, plot2, plot3, plot4, ncol=2, nrow=2)
  dev.off()
  
 
  # PID/ Table 1/ Figure 2
  # ## Table ###
  # out1<-rbind(coef(summary(m1)) [c(2:6),1:2], c(NA,NA), c(NA,NA), c(NA,NA), c(NA,NA)) #ANES
  # 
  # 
  # out2<-rbind(coef(summary(m2)) [c(2:10),1:2]) #ANES
  # results1<-as.matrix(
  #   cbind(out1, out2
  #   ))
  # 
  # row.names(results1)<-c("Operational Ideology", "In Party Negative", "Out Party Negative", "In Party Positive", "Out Party Positive", "Out Party Negative x Ideology", "In Party Negative x Ideology", "Out Party Positive x Ideology", "In Party Positive x Ideology")
  # require(xtable)                  
  # xtable(results1, digits=3,
  #        caption="Ideological consistency (Left column) and ideological alignment (right column). Left column is logit; right is ordered logit. Entries are point estimates and standard errors. 
  #        Entries in bold indicate a slope that is two times the estimated standard error. Ideology is operational ideology in the alignment column.
  #        All models control for religion, racial resentment, gender, education, and age. The full models
  #        may be found in the appendix.")
  # 
  #################################################################################
  #################################################################################
  # Affective Polarization ##
  #################################################################################
  #################################################################################
  davis.data$pid.3<-recode(davis.data$pid, "1='Democrat'; 2='Democrat'; 3='Democrat'; 4='Independent'; 5='Republican'; 6='Republican'; 7='Republican'") #Lean
  davis.data<-subset(davis.data, !is.na(pid.3))  ### Ignore the non-respondents for pid
  
  davis.data$anger.out<-NA
  davis.data$anger.out[davis.data$pid.3=="Democrat"]<-davis.data$rep.anger[davis.data$pid.3=="Democrat"]
  davis.data$anger.out[davis.data$pid.3=="Republican"]<-davis.data$dem.anger[davis.data$pid.3=="Republican"]
  davis.data$anger.out[davis.data$pid.3=="Independent"]<-rowMeans(cbind(davis.data$dem.anger[davis.data$pid.3=="Independent"],data$rep.anger[data$pid.3=="Independent"]), na.rm=T)

  
  davis.data$fear.out<-NA
  davis.data$fear.out[davis.data$pid.3=="Democrat"]<-davis.data$rep.afraid[davis.data$pid.3=="Democrat"]
  davis.data$fear.out[davis.data$pid.3=="Republican"]<-davis.data$dem.afraid[davis.data$pid.3=="Republican"]
  davis.data$fear.out[davis.data$pid.3=="Independent"]<-rowMeans(cbind(davis.data$dem.afraid[davis.data$pid.3=="Independent"],data$rep.afraid[data$pid.3=="Independent"]), na.rm=T)
    
  davis.data$disgust.out<-NA
  davis.data$disgust.out[davis.data$pid.3=="Democrat"]<-davis.data$rep.disgust[davis.data$pid.3=="Democrat"]
  davis.data$disgust.out[davis.data$pid.3=="Republican"]<-davis.data$dem.disgust[davis.data$pid.3=="Republican"]
  davis.data$disgust.out[davis.data$pid.3=="Independent"]<-rowMeans(cbind(davis.data$dem.disgust[davis.data$pid.3=="Independent"],data$rep.disgust[data$pid.3=="Independent"]), na.rm=T)
  
  
  davis.data$worried.out<-NA
  davis.data$worried.out[davis.data$pid.3=="Democrat"]<-davis.data$rep.worried[davis.data$pid.3=="Democrat"]
  davis.data$worried.out[davis.data$pid.3=="Republican"]<-davis.data$dem.worried[davis.data$pid.3=="Republican"]
  davis.data$worried.out[davis.data$pid.3=="Independent"]<-rowMeans(cbind(davis.data$dem.worried[davis.data$pid.3=="Independent"],data$rep.worried[data$pid.3=="Independent"]), na.rm=T)
  
  davis.data$worried.in<-NA
  davis.data$worried.in[davis.data$pid.3=="Democrat"]<-davis.data$dem.worried[davis.data$pid.3=="Democrat"]
  davis.data$worried.in[davis.data$pid.3=="Republican"]<-davis.data$rep.worried[davis.data$pid.3=="Republican"]
  davis.data$worried.in[davis.data$pid.3=="Independent"]<-rowMeans(cbind(davis.data$dem.worried[davis.data$pid.3=="Independent"],data$rep.worried[data$pid.3=="Independent"]), na.rm=T)
  
  
  davis.data$anger.in<-NA
  davis.data$anger.in[davis.data$pid.3=="Democrat"]<-davis.data$dem.anger[davis.data$pid.3=="Democrat"]
  davis.data$anger.in[davis.data$pid.3=="Republican"]<-davis.data$rep.anger[davis.data$pid.3=="Republican"]
  davis.data$anger.in[davis.data$pid.3=="Independent"]<-rowMeans(cbind(davis.data$dem.anger[davis.data$pid.3=="Independent"],data$rep.anger[data$pid.3=="Independent"]), na.rm=T)
  
  davis.data$fear.in<-NA
  davis.data$fear.in[davis.data$pid.3=="Democrat"]<-davis.data$dem.afraid[davis.data$pid.3=="Democrat"]
  davis.data$fear.in[davis.data$pid.3=="Republican"]<-davis.data$rep.afraid[davis.data$pid.3=="Republican"]
  davis.data$fear.in[davis.data$pid.3=="Independent"]<-rowMeans(cbind(davis.data$dem.afraid[davis.data$pid.3=="Independent"],data$rep.afraid[data$pid.3=="Independent"]), na.rm=T)
  
  davis.data$disgust.in<-NA
  davis.data$disgust.in[davis.data$pid.3=="Democrat"]<-davis.data$dem.disgust[davis.data$pid.3=="Democrat"]
  davis.data$disgust.in[davis.data$pid.3=="Republican"]<-davis.data$rep.disgust[davis.data$pid.3=="Republican"]
  davis.data$disgust.in[davis.data$pid.3=="Independent"]<-rowMeans(cbind(davis.data$dem.disgust[davis.data$pid.3=="Independent"],data$rep.disgust[data$pid.3=="Independent"]), na.rm=T)
  
  davis.data$excited.out<-NA
  davis.data$excited.out[davis.data$pid.3=="Democrat"]<-davis.data$rep.excited[davis.data$pid.3=="Democrat"]
  davis.data$excited.out[davis.data$pid.3=="Republican"]<-davis.data$dem.excited[davis.data$pid.3=="Republican"]
  davis.data$excited.out[davis.data$pid.3=="Independent"]<-rowMeans(cbind(davis.data$dem.excited[davis.data$pid.3=="Independent"],data$rep.excited[data$pid.3=="Independent"]), na.rm=T)
  
  davis.data$hopeful.out<-NA
  davis.data$hopeful.out[davis.data$pid.3=="Democrat"]<-davis.data$rep.hopeful[davis.data$pid.3=="Democrat"]
  davis.data$hopeful.out[davis.data$pid.3=="Republican"]<-davis.data$dem.hopeful[davis.data$pid.3=="Republican"]
  davis.data$hopeful.out[davis.data$pid.3=="Independent"]<-rowMeans(cbind(davis.data$dem.hopeful[davis.data$pid.3=="Independent"],data$rep.hopeful[data$pid.3=="Independent"]), na.rm=T)
  
  davis.data$excited.in<-NA
  davis.data$excited.in[davis.data$pid.3=="Democrat"]<-davis.data$dem.excited[davis.data$pid.3=="Democrat"]
  davis.data$excited.in[davis.data$pid.3=="Republican"]<-davis.data$rep.excited[davis.data$pid.3=="Republican"]
  davis.data$excited.in[davis.data$pid.3=="Independent"]<-rowMeans(cbind(davis.data$dem.excited[davis.data$pid.3=="Independent"],data$rep.excited[data$pid.3=="Independent"]), na.rm=T)
  
  davis.data$hopeful.in<-NA
  davis.data$hopeful.in[davis.data$pid.3=="Democrat"]<-davis.data$dem.hopeful[davis.data$pid.3=="Democrat"]
  davis.data$hopeful.in[davis.data$pid.3=="Republican"]<-davis.data$rep.hopeful[davis.data$pid.3=="Republican"]
  davis.data$hopeful.in[davis.data$pid.3=="Independent"]<-rowMeans(cbind(davis.data$dem.hopeful[davis.data$pid.3=="Independent"],data$rep.hopeful[data$pid.3=="Independent"]), na.rm=T)
  
  davis.data$anxiety.out<-rowMeans(cbind(davis.data$anger.out, davis.data$fear.out, davis.data$worried.out, davis.data$digusted.in), na.rm=T) 
  davis.data$anxiety.in<-rowMeans(cbind(davis.data$anger.in, davis.data$fear.in, davis.data$worried.in, davis.data$digusted.in), na.rm=T) 
  davis.data$enthusiasm.in<-rowMeans(cbind(davis.data$excited.in, davis.data$hopeful.in), na.rm=T)
  davis.data$enthusiasm.out<-rowMeans(cbind(davis.data$excited.out, davis.data$hopeful.out), na.rm=T)
  
  davis.data$intAnxIn<-davis.data$ideology*davis.data$anxiety.in
  davis.data$intAnxOut<-davis.data$ideology*davis.data$anxiety.out
  davis.data$intEnthIn<-davis.data$ideology*davis.data$enthusiasm.in
  davis.data$intEnthOut<-davis.data$ideology*davis.data$enthusiasm.out
  
  
  ### Code Affective Polarization ###
  
  davis.data$marry.out<-NA
  davis.data$marry.out[davis.data$pid.3=="Democrat"]<-davis.data$marry.rep[davis.data$pid.3=="Democrat"]
  davis.data$marry.out[davis.data$pid.3=="Republican"]<-davis.data$marry.dem[davis.data$pid.3=="Republican"]
  
  davis.data$marry.in<-NA
  davis.data$marry.in[davis.data$pid.3=="Democrat"]<-davis.data$marry.dem[davis.data$pid.3=="Democrat"]
  davis.data$marry.in[davis.data$pid.3=="Republican"]<-davis.data$marry.rep[davis.data$pid.3=="Republican"]
  
  davis.data$move.out<-NA
  davis.data$move.out[davis.data$pid.3=="Democrat"]<-davis.data$move.rep[davis.data$pid.3=="Democrat"]
  davis.data$move.out[davis.data$pid.3=="Republican"]<-davis.data$move.dem[davis.data$pid.3=="Republican"]
  
  davis.data$move.in<-NA
  davis.data$move.in[davis.data$pid.3=="Democrat"]<-davis.data$move.dem[davis.data$pid.3=="Democrat"]
  davis.data$move.in[davis.data$pid.3=="Republican"]<-davis.data$move.rep[davis.data$pid.3=="Republican"]
  
  
  
  ## move.neighbor
  
 summary(polr(as.factor(marry.out)~
           ideology+anxiety.in+anxiety.out+enthusiasm.in+enthusiasm.out+
           catholic+no.religion+other.religion+
           female+education+news+non.white+
           age, data=davis.data))
 summary(polr(as.factor(move.out)~
                ideology+anxiety.in+anxiety.out+enthusiasm.in+enthusiasm.out+
                catholic+no.religion+other.religion+
                female+education+news+non.white+
                age, data=davis.data))
  
 
 summary(polr(as.factor(move.neighbor)~
                ideology+anxiety.in+anxiety.out+enthusiasm.in+enthusiasm.out+
                catholic+no.religion+other.religion+
                female+education+news+non.white+
                age, data=davis.data))
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
#### Plot Marginal Effects ######  
  
model1<-lm(policy.conservatism~ideology3+anxiety+enthusiasm+intAnx+intEnth+
             catholic+no.religion+other.religion+racial.resentment+
             north.central+northeast+west+female+education+
             age+cross.section+news, data=d1)
marginal.ols<-function(output, cov, simulations){
    require(MASS)
    sim.beta <- mvrnorm(simulations,coef(output),vcov(output)) ##Simulate Draws from a multivariate normal
    #design<-c(0,0,1,1)  #Pass the correct design matrix
    c.interval<-sim.beta
    conditional.effect<-c.interval[,2]+c.interval[,5]*cov
    lower.bound<-quantile(conditional.effect, probs=0.025)
    upper.bound<-quantile(conditional.effect, probs=.975)
    t<-cbind(mean(conditional.effect), lower.bound, upper.bound )
    return(t)  ##Return the marginal effect for the moderated variable
  }
  ##Marginal effect for ideology
x.range<-c(0:100)/100
  tt<-cbind(NA, NA, NA)
  for(i in 1:length(x.range)){
    tt<-rbind(tt, marginal.ols(model1, x.range[i], 10000 ))
  }
tt<-tt[-1,]
plot.data<-as.data.frame(cbind(tt, x.range))
names(plot.data)<-c("mean", "min", "max", "Fear")
  ##Marginal effect for policy conservatism
  
  plot<-
    ggplot(data=plot.data) + 
    geom_ribbon(aes(x=Fear, ymin=min, ymax=max), fill="grey", alpha=0.5)+
    geom_line(aes(x=Fear, y=mean))+
    theme_bw() +
    theme(panel.background=element_rect(fill="#F0F0F0")) +
    theme(plot.background=element_rect(fill="#F0F0F0")) +
    theme(panel.border=element_rect(colour="#F0F0F0")) +
    # Format the grid
    theme(panel.grid.major=element_line(colour="#D0D0D0",size=.25)) +
    (axis.ticks=element_blank())+
    ggtitle("Marginal Effect of Symbolic Ideology on Operational Ideology") +
    theme(plot.title=element_text(face="bold",hjust=-.08,vjust=2,colour="#3C3C3C",size=13)) +
    theme(axis.text.x=element_text(size=11,colour="#535353",face="bold")) +
    theme(axis.text.y=element_text(size=11,colour="#535353",face="bold")) +
    theme(axis.title.y=element_text(size=11,colour="#535353",face="bold",vjust=1.5)) +
    theme(axis.title.x=element_text(size=11,colour="#535353",face="bold",vjust=-.5)) +
    scale_y_continuous("Symbolic Ideology", limits=c(0,4))+
    scale_x_continuous("Fear")
  
  
  pdf("/users/chrisweber/Dropbox/Project Folder/projects/Nick Davis_Polarization Emotion/paper/graph2.pdf", width=8, height=7)
plot
  dev.off()
  
##########################################################################
### Two things need to be done here:

pdf("/users/chrisweber/Dropbox/Project Folder/projects/Nick Davis_Polarization Emotion/paper/graph3.pdf", width=8, height=7)
  hist(d1$policy.conservatism, prob=T, xlab="Operational Ideology", main="Operational Ideology (by Symbolic Ideology)", ylim=c(0,0.6))
  lines(density(d1$policy.conservatism[which(d1$id.3cat=="Conservative")], adjust=1.5), ylim=c(0,1), ylab="Density", xlab="Operational Ideology", main="", col="red")
  lines(density(d1$policy.conservatism[which(d1$id.3cat=="Moderate")], adjust=1.5), col="purple")
  lines(density(d1$policy.conservatism[which(d1$id.3cat=="Liberal")], adjust=1.5), col="blue")
  legend(-2.1,0.6, legend=c("Symbolic Conservative", "Symbolic Moderate", "Symbolic Liberal"), lty=1,  col=c("red", "purple", "blue"), cex=0.7)
dev.off()
  
range(d1$policy.conservatism, na.rm=T)
#### Fill out 3x3 table #####

d1$intFear.moderate<-d1$moderate.s*d1$fear
d1$intFear.conservative<-d1$conservative.s*d1$fear
  
summary(polr(as.factor(group)~moderate.s+conservative.s+anger+
               fear+hope+pride+intFear.moderate+
               intFear.conservative+catholic+no.religion+
               other.religion+racial.resentment+north.central+
               northeast+west+female+education+age,data=d1))

model1<-polr(as.factor(group)~moderate.s+conservative.s+anger+
               fear+hope+pride+intFear.moderate+
               intFear.conservative+catholic+no.religion+
               other.religion+racial.resentment+north.central+
               northeast+west+female+education+age,data=d1)


pred.function.table<-function(output, moderate=1, conservative=1, fear=0, group="Conservative"){
  t<-model.matrix(output)[,-1]  # The data
  model.matrix<-cbind(
    ifelse(moderate==1, 1, 0),
    ifelse(conservative==1, 1, 0),
    mean(t[,3], na.rm=T), # anger
    fear,  #Fear
    mean(t[,5], na.rm=T), # hope
    mean(t[,6], na.rm=T), # pride
    ifelse(moderate==1, 1, 0)*fear,
    ifelse(conservative==1, 1, 0)*fear,
    min(t[,9], na.rm=T), #Catholic
    min(t[,10], na.rm=T), #None
    min(t[,11], na.rm=T), #Other
    mean(t[,12], na.rm=T), # Racial Resentment
    min(t[,13], na.rm=T), # North Centra
    min(t[,14], na.rm=T), # Northeast
    min(t[,15], na.rm=T), # West
    min(t[,16], na.rm=T), # Female
    max(t[,17], na.rm=T), # Education
    mean(t[,18], na.rm=T) #Age
  )
  beta.sim<-mvrnorm(1000, c(output$coefficients, output$zeta), vcov(output)) ##Draw samples from multivariate distrbution
  g1<-cbind(model.matrix, -1, 0)%*%t(beta.sim)
  g2<-cbind(model.matrix, 0, -1)%*%t(beta.sim)
  pr.1<-1-plogis(g1)
  pr.2<-plogis(g1)-plogis(g2)
  pr.3<-plogis(g2)
  if(group=="Moderate"){
    cat("The min, mean, max Mod prediction is", c(quantile(pr.2, probs=0.025), quantile(pr.2, probs=0.5), quantile(pr.2, probs=0.975)))
  }
  if(group=="Liberal"){
    cat("The min, mean, max Lib prediction is", c(quantile(pr.1, probs=0.025), quantile(pr.1, probs=0.5), quantile(pr.1, probs=0.975)))
  }
  if(group=="Conservative"){
    cat("The min, mean, max Con prediction is", c(quantile(pr.3, probs=0.025), quantile(pr.3, probs=0.5), quantile(pr.3, probs=0.975)))
  }
}  

pred.function.table(model1, moderate=0, conservative=0, fear=0, "Liberal")
pred.function.table(model1, moderate=0, conservative=0, fear=0, "Moderate")
pred.function.table(model1, moderate=0, conservative=0, fear=0, "Conservative")

pred.function.table(model1, moderate=1, conservative=0, fear=0, "Liberal")
pred.function.table(model1, moderate=1, conservative=0, fear=0, "Moderate")
pred.function.table(model1, moderate=1, conservative=0, fear=0, "Conservative")

pred.function.table(model1, moderate=0, conservative=1, fear=0, "Liberal")
pred.function.table(model1, moderate=0, conservative=1, fear=0, "Moderate")
pred.function.table(model1, moderate=0, conservative=1, fear=0, "Conservative")

pred.function.table(model1, moderate=0, conservative=0, fear=1, "Liberal")
pred.function.table(model1, moderate=0, conservative=0, fear=1, "Moderate")
pred.function.table(model1, moderate=0, conservative=0, fear=1, "Conservative")

pred.function.table(model1, moderate=1, conservative=0, fear=1, "Liberal")
pred.function.table(model1, moderate=1, conservative=0, fear=1, "Moderate")
pred.function.table(model1, moderate=1, conservative=0, fear=1, "Conservative")

pred.function.table(model1, moderate=0, conservative=1, fear=1, "Liberal")
pred.function.table(model1, moderate=0, conservative=1, fear=1, "Moderate")
pred.function.table(model1, moderate=0, conservative=1, fear=1, "Conservative")

range(d1$fear, na.rm=T)

pred.function.table(model1, moderate=0, conservative=0, fear=1, "Conservative")
